%% Output statistic generation function with multidimensional random variables as input
%% ============================ [function description] ============================
%% Prerequisite files.
%polCoeCal.m
%Inputs.
%inputData------ input values for matching points, matrix
%outputData----- output value of the matching point, matrix
%order---------- polynomial order, value

%MCnumber------- number of Monte Carlo simulation points, numeric
%output.
%PCmean--------- mean value by orthogonal polynomial property, numerical
%PCvar---------- variance by orthogonal polynomial property, numerical
%MCmean--------- mean by Monte Carlo and polynomial, numerical
%MCvar---------- variance by Monte Carlo and polynomial, values
%MCdata--------- Monte Carlo method output points, matrix
%MCinput-------- Monte Carlo method input points, matrix
%% ============================ [function body] ============================
function [PCmean,PCvar,MCmean,MCvar,MCdata,MCinput]=GetStatistics(inputData,outputData,order,MCnumber)
[c,y]=polCoeCal(inputData,outputData,order);%Get Ermitian polynomial coefficients, polynomial expressions
Ds=length(inputData(1,:));
MCinput=randn(MCnumber,Ds);
for i=1:Ds
    y=subs(y,sym(['x',num2str(i)]),MCinput(:,i));%Bring in discrete normal random variables
end
MCdata=double(vpa(y,8)); %Numerical symbolic variables
MCmean=mean(MCdata);     % Monte Carlo mean
MCvar=var(MCdata);       %Monte Carlo variance
PCmean=c(1);             %mean by orthogonal polynomial, with the first term of the polynomial being the mean of the random variable
PCvar=0;
for i=2:length(c)       %-|
    PCvar=PCvar+c(i)^2; %-|---Variance by orthogonal polynomial, order of all terms except first term
end                     %-| the sum of the squares of the products of the factorials and coefficients of